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Abstract 



We investigate, using mean-field theory and simulation, the ef- 
fect of asymmetry on the critical behavior and probability density 
of Bak-Sneppen models. Two kinds of anisotropy are investigated: 
(i) different numbers of sites to the left and right of the central 
(minimum) site are updated and (ii) sites to the left and right of 
the central site are renewed in different ways. Of particular in- 
terest is the crossover from symmetric to asymmetric scaling for 
weakly asymmetric dynamics, and the collapse of data with differ- 
ent numbers of updated sites but the same degree of asymmetry. 
All non-symmetric rules studied fall, independent of the degree of 
asymmetry, in the same universality class. Conversely, symmet- 
ric variants reproduce the exponents of the original model. Our 
results confirm the existence of two symmetry-based universality 
classes for extremal dynamics. 
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I. INTRODUCTION 



Nature exhibits scale invariance in a variety of settings. Such behavior is 
often associated with power law distributions of the phenomenon of interest, for 
example earthquake sizes [1,2], rain intensities and drought durations [3-5] and 
physiological and morphological quantities [6,7]. In recent decades, physicists 
have sought the physical origin of such laws [2,5,7]. A concept introduced 
to partially explain the ubiquity of scale-invariant phenomena in nature is 
so-called self-organized criticality or SOC [2,8]. 

The Bak-Sneppen (BS) model [8,9] was proposed to explain mass extinc- 
tions observed in the fossil record, and has attracted much attention as a 
prototype of SOC under extremal dynamics [10]. The model has been studied 
through various approaches, including simulation [11-14], theoretical analy- 
sis [15-17], probabilistic analysis (run time statistics) [18,19], renormalization 
group [20,21], field theory [22] and mean-field theory [9,10,23-27]; it was re- 
cently adapted to model experimental data on bacterial populations [28,29]. 
Applications of the model in evolution studies are reviewed in [30] . Asymmet- 
ric versions of the model were studied in Refs. [31] and [32]. In this paper, we 
study how varying the degree of symmetry affects the stationary distribution 
and the critical behavior, using mean-field theory and numerical simulation. 

In the evolutionary interpretation of the BS model, each site i represents a 
biological species, and bears a real- valued variable Xi representing its "fitness". 
The larger Xi, the better adapted this species is to its environment and so the 
more likely it is to survive. At each time step, the site with the smallest Xi 
and its nearest neighbors are replaced with randomly chosen values. The re- 
placement of Xi represents extinction of the less-fit species and the appearance 
of a new one, while the substitution of the neighboring variables with new 
random values may be interpreted as a sudden unpredictable change in fitness 
when an interdependent species goes extinct and a new species colonizes its 
niche. Selection, at each step, of the global minimum of the {xj} ("extremal 
dynamics" ) represents a highly nonlocal process, and would appear to require 
an external agent with complete information regarding the state of the system 
at each moment. 

Due to the extremal dynamics, this system exhibits scale-invariance in the 
stationary state, in which several quantities display power-law behavior [8]. 
Simulations show that the stationary distribution of barriers follows a step 
function, being zero (in the infinite-size limit) for x < ~ 0.66702(3) [12]. 
Relaxing the extremal condition leads to a smooth probability density and loss 
of scale invariance [10,26,27]. Datta et. al have also shown that the scaling 
behavior is sensitive to the number of sites updated at each step (i.e., updating 
only the minimum site, or the minimum and next-to-minimum as well) [33]. 

In this paper we investigate the effect of symmetry in the updating rule. 
Section II introduces the models, which are then analyzed using mean-field 
like approaches in Sec. III. In Sec. IV we present simulation results; our 
conclusions are summarized in Sec. V. 
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II. MODELS 



The Bak-Sneppen model [8] is defined as follows. Consider a dimensional 
lattice with sites and periodic boundaries. In the evolutionary interpre- 
tation of the model, each site i represents a species, and bears a real-valued 
variable Xi{t) representing its "fitness" to survive, so that Xi{t) may be termed 
a "barrier to extinction". The initial values of the barriers are independently 
drawn from a uniform distribution on the interval [0,1). At time 1, the site 
m bearing the minimum of all the numbers {xiiQ)} is identified, and it, along 
with its 2d nearest neighbors, are given new random values, again drawn in- 
dependently from the interval [0,1). (In the one dimensional case considered 
here this amounts to: XmiX) = fj, Xra+i{^) = f]' ■, and Xm-iiX) = f]" ■, where 77, 
77', and Tj" are independent and uniformly distributed on [0,1); for \j — m| > 1, 
Xj{l) — Xj{^).) At step 2 this process is repeated, with m representing the site 
with the global minimum of the variables {xj(l)}, and so on. 

We now define several one-dimensional variants of the BS model, which 
differ from the original in the number and/or position of neighbors which are 
updated at each time step, or in the way that the barriers Xi evolve. In the 
'generalized' or 'BSab' variant, we replace the site m bearing the minimum of 
the {xi\ plus a neighbors on the left side and h neighbors on the right with 
independent random numbers. If a = 6 = 1 we recover the original model 
(BSll); if a = and 6 = 1 we have the anisotropic BS model (BSOl) studied 
in [31,32]; if a 7^ 6 we obtain modified BS models with asymmetric dynamics. 
These BSab variants of the Bak-Sneppen model were also studied in [31]. In the 
second variant, the site M bearing the maximum of the {xj} and its two nearest 
neighbors are updated according to the rules: XM^t + 1) = XM'it + 1) — r]' 
and XM"{t + 1) = [xM"{t)f- Here M' = M + a and M" = M - a, where a 
is +1 with probability p and —1 otherwise. We shall refer to this variant as 
the peripheral square model with variable anisotropy. (For p = 1/2 this is the 
peripheral square model studied in [27].) 

The motivation for studying these variants is threefold. First it is of interest 
to examine the effect of (i) the symmetry of the dynamics and (ii) replacement 
of barriers with a deterministic function /(,t) instead of random numbers, 
on the critical behavior of the model. Secondly, we study corrections to the 
power-laws, such as finite-size effects, and find that these corrections are large 
for small asymmetries in agreement with [31]. Finally, since the precise form 
of the dynamics in a specific setting (e.g., evolution) is generally unknown, 
and probably is quite different from that of the original model, it is of interest 
to test the robustness of the results reported for the original model. It is even 
possible that some of the variants considered approximate a given evolutionary 
process more closely than the original. In particular, if certain pairs of species 
[i and j, say) stand in a predator-prey relationship, one would not expect the 
extinction of i to have the same effect on j as the extinction of j has on i; in 
this situation an asymmetric interaction appears more reasonable. 
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III. MEAN-FIELD THEORY 



We develop a mean-field theory for the BSab variants, along the lines of 
Refs. [10,26]. To begin, we introduce a flipping rate of Te'^^^ at site z, where 
is a characteristic time, irrelevant to stationary properties, and which we 
set equal to one (F = 1). This regularized system is the 'finite-temperature' 
model [10,26,27], in which all sites have a nonzero probability of being up- 
dated at any time, in contrast to the extremal dynamics of the original model. 
The extremal condition is recovered in the zero-temperature limit, /3 — > oo; a 
regularized flipping rate facilitates analysis of the model. 

The probabihty density p{x, t) satisfies: 



dt 



= -e-^M^, t)-±[' e->,(x, y, t)dy -j^C e-^^p,(x, y, t)dy 

+ n [\-^^p{y,t)dy, (1) 
Jo 



where n = a + 6-|-lis the number of sites that are updated at each step, 
Pj{x, y, t) is the joint density for sites and j and p{y, t) is the one-site marginal 
density. (We assume translation and reflection invariance, which is expected to 
hold at any time, if the initial distribution possesses these properties.) Invoking 
the mean- field factorization Pj{x, y, t) — p{x, t)p{y, t), we find: 

^P^ = -p{x,t)[e-(''' + in-l)im+nI{P) , (2) 

where 

m= f\-^yp{y,t)dy (3) 
represents the overall flipping rate. In the stationary state we have 

= {n - 1)1+ e-^^ ■ 

Multiplying by e~^^ and integrating over the range of x, we find — 
(e(»-i)/3/» - i)/[(n - 1)6^^(1 - e-^/")] and thus 

- _ 1) 1 _ e-{n-iWn + e-P^ie^/^ - 1) ' ^ ^ 

In the limit /3 — > oo this solution becomes a step function: 

77 

Pst{x) = -^^e{x-l/n)eil-x) . (6) 

Thus, the mean-field approach predicts a step-function singularity for the prob- 
ability density with the critical barrier at x* — 1/n. This result is independent 
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of the symmetry of the dynamics, and we conclude that the mean-field approx- 
imation at the level of the one-site marginal density is insensitive to differences 
in symmetry. 

Define the anisotropy coefficient 

K = (7) 

n 

(Note that ka = for symmetric dynamics.) The mean-field threshold x* — 
1/n can be written as 

X* = (^^^^(1 - ^a) for b>a. (8) 

For fixed a, mean-field theory predicts that the threshold x* varies hnearly 
with the anisotropy coefficient ka- 



IV. SIMULATIONS 
A. Threshold values 

We study the generalized model for various values of a and b, corresponding 
to ka in the range 0—0.857. We estimate the probability density p{x) on the 
basis of a histogram of barrier frequencies, dividing [0,1] into 100 subintervals. 
After the system (a ring of = 2000 sites) relaxes to the stationary state-^, 
histograms are accumulated at intervals of time steps until a total of 10^ 
steps are performed. 

The threshold x*pf for a finite system can be calculated as follows [34] . In 
the limit A^ oo, the probability density p{x) is a step function, with p{x) = 
for X < X* and p{x) = C for x > x* ^ where C is a constant. Normalization 
implies that C{l — x*) = 1. This suggests that = 1 — 1/Cn can be regarded 
as the threshold of a finite system. In ref. [34] we performed simulations for 
various systems sizes and found that = A;A^-V^ with x*^ = 0.6672(2) 

and u = 1.40(1) in the original model (or BSll) and x^ = 0.7240(1) and 
u = 1.58(1) in the anisotropic BS model (or BSOl). Typically x*j^ — xl^ ^ 0.006 
for A^ = 2000 and, therefore, the threshold a;J^^2ooo sufficiently accurate for 
our present purpose, which is to investigate the effect of anisotropy on the 
threshold values. 

Simulation confirm that the stationary probability density is a step func- 
tion, in agreement with MFT. Fig. 1 shows p{x) for all variants having n = 5. 
The threshold values, listed in Table I, are however always larger than that 
predicted by MFT, x* — 1/n. MFT is exact for the random-neighbor version. 



^Typically, the system can be considered to be in the stationary state after !=a 10^ 
time steps. 



in which all sites arc considered neighbors, and in which (n — 1) randomly 
selected sites are updated in addition to the site with the minimum Xj. There 
is better agreement between the MFT prediction for x* and the simulation 
result for larger anisotropics (with n constant), and for larger n (with fixed 
ka). MFT predicts that x* — x*{ri), whereas in fact, for fixed n, the threshold 
is smaller, the larger the anisotropy. This may be understood in general terms 
as follows. With n fixed, the larger |a — 6|, the greater is the distance between 
updated sites and the minimum site. This makes the dynamics resemble the 
mean-field (random neighbor) case more closely, causing the threshold to tend 
toward the MFT value 1/n. Increasing n (with ka fixed) should have a similar 
effect. 

Figure 2 shows that if we fix a and vary b, the threshold decreases linearly 
with ka, as predicted by MFT, except for a = 0. The actual threshold values 
are, however, quite different from MFT results. For instance for a = 1, MFT 
predicts x* — (1/3) (1 — /;;„), while a fit to the simulation data yields x* — 
0.664(3) - 0.668(7)A;„. 

B. Critical exponents 

Several quantities display power-law behavior in the Bak-Sneppen model, 
and can be used to characterize the associated universality class [8,11,31]. In 
particular, we study the probability Pj{r) that successive updated sites are 
separated by a distance r. In the original model, one finds [31] 

Pj(r) ~ r--^, with 71 = 713 = 3.23(2) (9) 

(figures in parentheses denote uncertainties). On the other hand, in the 
anisotropic BS model (BSOl), vr = tt^ = 2.401(2) [32]. 

Based on simulations of the generalized models using 7 x 10^ time steps 
on lattices of 32000 sites we find that all variants with asymmetric dynamics 
belong to the universahty class of the anisotropic BS model. Fig. 3 displays 
Pj{t) for variants with different anisotropics. The corrections to the power 
laws are large for small asymmetries and therefore the asymptotic behavior is 
not observed in small lattices in these cases, because large values of r are not 
reached. 

We find that Pj{r) can be fitted with an expression representing a crossover 
between symmetric scaling at small r and asymmetric scaling at large r, 

Pj{r-) = A[r-^^ + {N -r)-^'']+ B[r-^s + {N -r)-^s] , (10) 

where A and B are constants and tt^ and tt^ are, respectively, the exponents of 
the anisotropic and the original models quoted above. The 'mirror' structure 
of Eq. (10) arises due to the periodic boundary conditions, which imply that 
Pj{'<') = Pj{N — f)- Eq. (10) provides a good fit to the simulation data 
for r > 10^ for all generalized models. Since i^s ^ '^a + 1, the expression 
Pj(r) = Ar-^^(1 + B/r) + A{N - r)-^^(l + B/{N - r)) also fits the data 
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well. Nevertheless, studies of weakly asymmetric models (e.g., BS(20)(21), for 
which ka = 0.024) indicate that the slope of Pj(r) (on log scales) approaches 
Tis-i as opposed to tta + 1, before the asymptotic behavior is attained. 

The functions Pj{r) for variants with the same anisotropy can be collapsed 
onto a master curve 

Pcoiiapse{ka,r') = n Pj{r/n) , (11) 

where r' — r/n. Figure 4 illustrates this collapse for five variants with ka — 1/6. 
A simple argument provides an intuitive understanding of this scaling function. 
Suppose site m is the extremal site at time t. Since renewed sites receive 
independent random numbers, all sites in the set Mit) = {m{t) — a,m{t) — 
a + 1, ...,m{t) + b} have the same probability of being the extremum at time 
t + 1, which implies (with a < b) 

2Pj{0) = Pj{l) = Pj{2) = ... = Pj{a), (12) 

so that Pj{r) exhibits a plateau for r = 1, a. 

The probability of event m{t + 1) G A/'(t), i.e., the extremal site belongs to 
the group of sites updated at the last time step, can be estimated as follows. 
We assume that, with a probability ^ 1, all sites outside of J\f{t) have Xi > x*. 
Then the probability that m{t + 1) G Af{t) is 1 - (1 - a;*)". In MFT we then 
have P{1) = ^[1 — {1 — X*)"] oc ^, consistent with the prefactor n in Eq. (11). 
(Observe that for the original BS model, 1 — (1 — x*)" ~ 0.96.) 

U b > a, then 

Pj{a + 1) ~ Pj(a + 2) ~ • • • ~ Pj{b) ~ ^Pj(l) (13) 

where the symbol '~' is used because there is a small probability of one of the 
distances r = a + 1, 6 being the extremum distance even if m(t + 1) ^ A/'(t). 
This explains the second plateau seen in Fig. 4. For fixed ka, a = [n{l — ka) — 
l]/2 and b = [n(l + ka) — l]/2 are essentially proportional to n, so that the 
rescaling of the argument in Eq. (11) collapses the plateaus. For r > b the 
probability Pj rapidly approaches a power law, which is also invariant under 
the rescaling of Eq. (11). Although the foregoing arguments are approximate, 
they provide an intuitive basis for the collapse seen in Fig. 4. 

Another quantity that obeys a power law in the Bak-Sneppen model is 
Pr{T), the probability that, in the stationary state, the global minimum site 
at time t was the extremal site most recently at time t — r. (r is the 'return 
time.') We simulate the BSab variants on a lattice of 16000 sites, using 7 x 10* 
time steps. The results are quite similar to those for the probability Pj{r). 
For weak anisotropics, the asymptotic behavior is observed only for large r, as 
shown in Fig. 5. In the limit r — > oo, Prir) oc t~^, where b — 1.40(1) for all 
anisotropic variants and b— 1.58(1) for all isotropic cases. 
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C. Peripheral square model with variable anisotropy 



The peripheral square model with variable anisotropy, defined in Sec. II, 
is a generalization of the model studied in [27], where the anisotropy now 
depends on a parameter p. We simulated this model (with N — 2000, 4000 
and 8000 sites) for p varying from the isotropic value {p = 0.5) to the maximally 
anisotropic one {p = 1.0). Our analysis shows that, in the limit L — oo, the 
stationary probability density consists of two regions (see Fig. 6): Pst{x) = 
for X > X* and Pst{x) 7^ for x < x*. (The rounding in the threshold region is 
simply a finite-size effect.) This behavior is in agreement with Ref. [27], where 
we found that renewing the barriers via x' = with a = 1/2 and 2, leads to a 
diversity of distributions with a discontinuity at the threshold. The divergence 
of Pst as a; in the present case can be understood on the basis of the mean- 
field theory developed in [27]. We conclude that the step-function distribution 
of the original BS model is not universal for self-organized criticality under 
extremal dynamics. In contrast, our results suggest that the step-singularity 
is universal. 

Figure 6 suggests that, in the limit a; 0, Pst{x) does not depend on the 
parameter p. Over restricted intervals Pst appears to follow a power law. For 
example, for x e [10~^,10~^], Pst is well described by a power law with an 
exponent that increases from 0.82 for p — 0.5 to 0.85 for p = 1.0. On the 
other hand, the mean-field solution, which appears to capture the qualitative 
behavior, does not follow a power law, diverging instead as J2'^=o 4~'^a;^ in 
the limit x — > 0. (We further note that although the MF solution accompanies 
the general trend of the data, it exhibits a series of step singularities in addition 
to the jump at x*. Several step discontinuities are in fact observed in Pst in 
simulations of the random-neighbor version [27], although not in the nearest- 
neighbor version.) 

The dependence of the threshold x* on p is shown in Fig. 7. Note that 
augmenting the anisotropy, x* increases almost linearly. Therefore, similarly 
to the BSab variant, the effect of increasing the anisotropy, while maintaining 
the number of updated sites constant, is to decrease the interval on which 
p{x) = 0. 

Figure 8 shows the probabihty density Pj{r) for various values of p. While 
the isotropic case {p — 0.5) preserves the exponent of the original BS model, 
all cases for p 7^ 0.5 exhibit the exponent of the anisotropic BS model. In 
the BSab variant studied above, asymmetry was due to different numbers of 
sites being renewed in each side of the extremal site. In the present case 
we have a different kind of asymmetry, namely, sites are updated in different 
ways on each side of the extremal one. We nevertheless find the same critical 
exponents as for the anisotropic BS model. This strengthens the evidence 
for the existence of two symmetry-based universality classes for models under 
extremal dynamics. 
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V. CONCLUSIONS 



We perform a detailed investigation of the effect of symmetry on the scahng 
behavior of the Bak-Sneppen model. All dynamics which preserve the reflec- 
tion symmetry of the original BS model possess the same critical exponents 
as the original model, while asymmetric dynamics lead to the exponents of 
the anisotropic BS model. Therefore, our work reinforces the evidence for two 
symmetry-based universality classes [31]. 

In order to obtain these results, we deflne modifled BS models, and study 
them via mean-field theory and simulation. In the generalized BS model (or 
BSab variant), the degree of asymmetry is quantified by the anisotropy coeffi- 
cient ka- Mean-field theory provides poor predictions for the threshold x*, but 
correctly predicts (i) that the stationary probability density is a step function, 
and (ii) that x* varies linearly with the number n of updated sites, as found in 
simulations on a one-dimensional lattice. The simulations also lead to several 
conclusions that go beyond MFT analysis: (i) the threshold value decreases 
as the degree of anisotropy is increased; (ii) for weak anisotropy we observe a 
crossover between symmetric and asymmetric scaling; (iii) for fixed anisotropy 
we find a scaling collapse of the probability Pj{r) that successive minimum 
sites have a separation r. 

In the peripheral square model with variable anisotropy, the degree of asym- 
metry is quantified by the parameter p that controls on which side of the ex- 
tremal site the neighbor is updated differently. Although this model is partially 
deterministic and has a different kind of asymmetry, we again encounter two 
symmetry-based universality classes. Moreover, we find that increasing the 
asymmetry, the region where p{x) = is reduced, as in the the BSab variant. 

These results, together with evidence for finite-size scaling [34], and con- 
nections with directed percolation [22], indicate that scaling in the BS model 
partakes of many of the characteristics associated with critical phenomena, 
both in and out of equilibrium. The particular distinguishing features of the 
Bak-Sneppen model and its variants appear to be associated with extremal 
dynamics. 
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TABLES 



TABLE I. Threshold (x^) for systems of A = 2000 sites for generahzed BSab 
models with various anisotropy coefficients (ka)- Figures in parentheses denote un- 
certainties. 



MODEL 


* 

^W=2000 




BSOl 


0.717(2) 


1/2 


BS02 


0.531(2) 


2/3 


BS03 


0.417(2) 


3/4 


BS04 


0.342(2) 


4/5 


BS05 


0.290(2) 


5/6 


BS06 


0.251(2) 


6/7 


BSll 


0.661(9) 





BS12 


0.501(3) 


1/4 


BS13 


0.398(2) 


2/5 


BS14 


0.329(3) 


1/2 


BS15 


0.280(2) 


4/7 


BS22 


0.466(9) 





BS23 


0.381(3) 


1/6 


BS24 


0.317(3) 


2/7 


BS25 


0.271(2) 


3/8 


BS26 


0.237(2) 


4/9 


BS27 


0.210(2) 


1/2 


BS33 


0.36(1) 
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FIGURES 
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FIG. 1. Stationary probability densities for the BSab with n = 5. 




FIG. 2. Threshold values as a function of the anisotropy coefficient ka- We fix a 
and change b to vary ka- Further explanation in text. 




r 



FIG. 3. Stationary probability Pj{r) for various BSab variants. 
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FIG. 4. Scaling collapse of Pj{r) for variants with ka = 1/6, see Eq. (11). 
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FIG. 5. Stationary probability Prir) for BSab variants. 
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FIG. 6. Stationary probability density p(a;) given by mean-field theory and sim- 
ulation (A^ = 8000 sites). The parameter p varies from the symmetric value (p=0.5) 
to the maximaly asymmetric one (p=1.0). 
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FIG. 7. Thresholds of the peripheral square models as a function of the param- 
eter p. 
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FIG. 8. Stationary probability -Pj(r) for the peripheral square model with dif- 
ferent values of p. Inset: Local slope vs. 1/r. (Note that the slope - exponent tt - 
approaches —tta = —2.401 in the limit r ^ oo in all cases with p ^ 0.5.) 
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